*==============================================================================*
*         Mean Pollution Emission Rates of New U.S. Vehicles, 1957-2019        *
*==============================================================================*


*---------------------------- setup 1957-1971 data ----------------------------*
u "replicationFiles/dataSTATA/combined/combined_aes73.dta", clear

ren emissions_new_CO_natl emissions_new_CO
ren emissions_new_HC_natl emissions_new_HC
ren emissions_new_NOX_natl emissions_new_NOX
keep emissions_new_CO emissions_new_HC emissions_new_NOX model_year

tempfile temp
sa "`temp.dta'", replace


u "replicationFiles/dataSTATA/combined/combined_newcars.dta", clear
preserve
	keep if model_year >= 2000 & model_year <= 2015
	replace sales_weight = round(sales_weight)
	collapse (mean) Wemissions_new_CO  = emissions_new_CO ///
					Wemissions_new_HC  = emissions_new_HC ///
					Wemissions_new_NOX = emissions_new_NOX [fw=sales_weight], ///
					by(model_year) fast
	tempfile wmean
	sa "`wmean.dta'", replace
restore

collapse (mean) emissions_new_CO emissions_new_HC emissions_new_NOX, by(model_year) fast

merge 1:1 model_year using "`wmean.dta'", assert(1 3) nogen

merge 1:1 model_year using "`temp.dta'", nogen 
sa "`temp.dta'", replace

foreach p in CO HC NOX {

	if inlist("`p'","CO","HC") loc xline 1967
	if "`p'" == "NOX" loc xline 1971
	
	if "`p'" == "CO"  loc ylab "85 30 10 3 1 0.4" 
	if "`p'" == "HC"  loc ylab "10 3 1 0.3 0.1 0.025"
	if "`p'" == "NOX" loc ylab "6 2 0.6 .2 .08 0.04 0.02"

	tw (connected emissions_new_`p' model_year, sort lcolor(blue) mcolor(blue) msymb(S)) ///
		(connected Wemissions_new_`p' model_year, sort lcolor(red) mcolor(red) msymb(Oh)), ///
		graphr(color(white)) ytit("Emission rate (grams per mile)") ///
		xtit("Model Year") yscale(noline) ///
		legend(order(1 "Unweighted mean" 2 "Weighted mean")) ///
		xlab(1957 1967 1977 1987 1997 2007 2020) ///
		ylab(`ylab') ///
		yscale(log) ///
		xline(`xline', lcolor(black)) 
		
	graph export "results/figures/f1_newcar_trend/`p'.wmf", replace
	graph export "$overleaf/figures/f1_newcar_trend/`p'.eps", replace

}

* for text summary of % change
foreach p in CO HC NOX {
	su emissions_new_`p' if model_year <= 1965
		loc rmean1 = r(mean)
	su emissions_new_`p' if model_year == 2020
		loc rmean2 = r(mean)		
	di `rmean2' / `rmean1'
}

*---------------------------------- CO2 ---------------------------------------*
* 1957-1974 from USEPA (1973)
* 1975-2020 from USEPA (2020)

insheet using "replicationFiles/dataRAW/mpg/1975_2020/table_export.csv", clear
ren *, l
keep if regulatoryclass == "All"
drop if modelyear == "Prelim. 2021"
assert vehicletype == "All"
assert real(productionshare) == 1
keep modelyear realworldmpg 
ren modelyear model_year
destring model_year, replace
ren realworldmpg mpg
su mpg if model_year == 1975
loc rmean_latedata = r(mean)
tempfile mpg_1975_2020
sa "`mpg_1975_2020.dta'", replace


import excel using "replicationFiles/dataRAW/mpg/1957_1974/fueleconomy_1957_1974.xlsx", clear firstrow
ren mpg_salesweighted_ftp72 mpg
su mpg if model_year == 1975
loc rmean_earlydata = r(mean)

* splice
replace mpg = mpg + (`rmean_latedata' - `rmean_earlydata')
drop if model_year == 1975

append using "`mpg_1975_2020.dta'"

g emissions_new_CO2 = 19.37 * 453 / mpg

tw (connected emissions_new_CO2 model_year, sort lcolor(blue) mcolor(blue) msymb(S)), ///
	graphr(color(white)) ytit("Emission rate (grams per mile)") ///
	xtit("Model Year") yscale(noline) ///
	xlab(1957 1967 1977 1987 1997 2007 2017) ///
	ylab(800 250 75 25 8 3 1) ///
	yscale(log) ///
	legend(order(1 "Weighted mean")) ///
	xline(1974.5, lcolor(black)) 
	graph export "results/figures/f1_newcar_trend/CO2.wmf", replace
	graph export "$overleaf/figures/f1_newcar_trend/CO2.eps", replace


su emissions_new_CO2 if model_year == 1957
	loc rmean1 = r(mean)
su emissions_new_CO2 if model_year == 2017
	loc rmean2 = r(mean)
	
di `rmean2' / `rmean1'
	
	